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Coherent  techniques  such  as  the  block  accumulating  coherent  integration  over  extended  interval  (BACIX)  have  recently  been 
proposed  to  extend  coherent  integration.  Although  efficient,  such  coherent  methods  may  still  be  too  expensive  except  for  high-end 
receivers  and  may  not  maintain  the  SNR  performance  when  there  are  large  frequency  changes  over  the  intended  integration 
interval. 

In  this  paper,  we  set  forth  a  novel  method  that  utilizes  the  semi-coherent  scheme  for  post-correlation  integration,  which  is  named  as 
“block  accumulating  semi-coherent  integration  of  correlations”  or  BASIC.  It  can  be  viewed  as  an  extension  of  the  BACIX 
algorithm.  Although  less  sensitive  than  coherent  integration,  semi-coherent  integration  based  on  inter-block  conjugate  products  is 
computationally  more  efficient.  In  addition,  it  can  handle  large  frequency  changes. 

The  BASIC  algorithm  is  first  formulated  in  the  paper.  Computer  simulation  results  are  then  presented  to  illustrate  operation  and 
performance  of  the  BASIC  algorithm  for  joint  estimation  of  the  initial  frequency,  chirping  rate  (frequency  rate),  bit  sync,  and  bit 
sign  pattern. 
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ABSTRACT 

This  paper  investigates  the  application  of  the  Wigner- 
Hough  transform  (WHT)  or  Wigner-Radon  transform 
(WRT)  and  their  variants  to  GPS  post-correlation 
integration.  For  an  acquisition  scheme  wherein  frequency 
search  step  is  500  Hz  and  code  search  step  is  a  half  chip, 
GPS  despreading  correlation  is  typically  carried  out  over 
1  millisecond.  Such  integration  at  the  rate  of  1  kHz  can 
tolerate  ±250  Hz  errors  in  frequency  and  ±14  chips  in  code 
phase  in  the  sense  of  3  dB  loss  of  correlation  power.  In 
many  applications,  it  is  desired  to  integrate  the  1  ms 
correlations  over  a  rather  long  period  of  time  so  as  to 
boost  the  signal  strength  while  averaging  out  noise  and 
interference.  When  the  integration  interval  extends  to  and 
beyond  20  ms,  the  signs  of  data  bits  have  to  be 
considered.  For  even  longer  intervals,  changes  in  signal 
frequency  due  to  motion-induced  Doppler  and/or  clock 
drift  become  important  and  have  to  be  taken  into  account. 
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Under  a  constant  acceleration,  the  complex  correlations 
are  subject  to  a  varying  Doppler  frequency,  which 
behaves  like  a  chirp  signal  (a  linear  frequency  modulation 
or  LFM).  The  Wigner-Ville  distribution  (WVD)  is  a  well 
known  method  to  estimate  instantaneous  frequency, 
which  appears  as  a  straight  line  in  the  time-frequency 
plane.  The  Hough  transform  (HT)  is  another  well-known 
method  to  detect  a  line,  which  can  be  used  to  integrate  the 
chirp  signal  power  along  the  line.  Similar  to  the  HT,  the 
Radon  transform  (RT)  is  another  formulation  that 
evaluates  integrals  along  straight  lines.  However,  the 
combination  of  the  WVD  and  HT  (called  the  Wigner- 
Hough  transform)  or  the  combination  of  the  WVD  and 
Radon  transform  (called  the  Wigner-Radon  transform), 
cannot  be  applied  directly  to  the  despread  correlations 
because  of  the  presence  of  unknown  data  bits. 

A  block- implemented  WHT/WRT  using  the  intra-block 
conjugate  products  is  set  forth  in  this  paper,  which  can  be 
used  to  estimate  the  initial  frequency,  chirping  rate 
(frequency  rate),  and  bit  sync.  Simulation  results  are 
presented  to  illustrate  the  operation  and  performance  of 
the  method. 

INTRODUCTION 

In  [Yang  et  al.,  2007],  six  post-correlation  integration 
schemes  are  comparatively  studied  with  computer 
simulation.  These  six  schemes  are  (1)  Ideal  Coherent,  (2) 
Practical  Coherent  with  FFT,  (3)  Non-Coherent,  (4)  Semi- 
Coherent  for  First  Lag  with  FFT,  (5)  Semi-Coherent  for 
First  Lag,  and  (6)  Semi-Coherent  up  to  First  N/2  Lags, 
which  are  illustrated  in  Figure  1. 


Ideal  Coherent 
Integration 


Practical  Coherent 
Integration  with  FFT 


Non-Coherent 

Integration 

Semi-Coherent 
Integration  for 
First  Lag  with  FFT 

Semi-Coherent 
Integration  for 
First  Lag 


Semi-Coherent 
Integration  for 
First  K  Lags 


Fig.  1  -  Block  Diagram  of  Six  Integration  Schemes 

When  the  change  in  carrier  frequency  is  small,  the  six 
integration  schemes  are  ranked  for  the  SNR  required  to 
achieve  the  same  probability  of  detection  for  a  given 
probability  of  false  alarm.  The  performance  ranking  has 
the  following  order: 


Ideal  Coherent  > 

Practical  Coherent  with  FFT  > 

Semi-Coherent  up  to  First  N/2  Lags  > 

Semi-Coherent  for  First  Lag  > 

Semi-Coherent  for  First  Lag  with  FFT  > 

Non-Coherent 

However,  the  performance  of  four  integration  schemes, 
except  for  the  coherent  and  non-coherent  schemes,  starts 
to  deteriorate  when  the  change  in  carrier  frequency  is 
large,  as  caused  by  motion-induced  Doppler  frequency 
shift  and  clock  instability  and  drift.  The  effect  is 
especially  pronounced  when  a  long  integration  interval  is 
used  in  order  to  cumulate  the  enough  signal  strength 
while  attenuating  noise. 

One  way  to  model  large  change  in  frequency  is  to  use  a 
quadratic  phase  function  or  a  linear  frequency  modulation 
(LFM)  where  the  frequency  change  rate  is  constant.  LFM 
is  widely  used  in  radar  signal  design,  which  produces  the 
so-called  chirp  signal  with  wide  instantaneous  bandwidth 
to  improve  range  resolution  [Wirth,  2001;  Yang  and 
Blasch,  2007].  It  can  also  be  used  to  model  a  variety  of 
phenomena  encountered  in  radar  signal  processing.  One 
example  is  the  phase  variation  in  the  cross  range  of  a 
synthetic  aperture  radar  (SAR).  It  is  fairly  accurate  to 
model  short-term  phase  change  with  a  chirp  signal. 

Since  the  chirping  rate  is  unknown  but  within  a  certain 
range,  one  technique  would  call  for  a  bank  of  processors, 
each  tuned  to  a  discretized  chirping  rate  in  the  range  (thus 
dechirping)  [Li,  1987].  Each  of  the  parallel  processors  can 
implement  one  of  the  six  integration  schemes  above 
mentioned.  This  is  tantamount  to  a  two-dimensional  (2D) 
search,  one  in  the  chirping  rate  and  the  other  in  the  initial 
frequency  with  FFT  used  for  integration.  The  use  of  an 
adaptive  fdter  to  implement  2D  search  is  reported  in 
[Wang,  Xia,  and  Chen,  1993]. 

There  are  many  methods  available  in  the  digital  signal 
processing  literature  for  estimating  a  chirp  signal’s 
parameters  ( i.e .,  the  initial  phase,  initial  frequency,  and 
chirping  rate)  and  one  method  of  particular  interest  is  the 
Wigner-Hough  transform  (WHT)  [Barbarossa,  1996].  In 
this  method,  the  Wigner-Ville  distribution  (WVD)  is  used 
to  represent  the  signal  energy  in  the  time-frequency  plane 
while  the  Hough  transform  (HT)  integrates  the  energy 
belonging  to  a  chirp  signal,  which  is  distributed  along  a 
line  in  the  time-frequency  plane.  Similarly,  the  Radon 
transform  can  be  used  in  place  of  the  Hough  transform  to 
integrate  the  chirp  energy  along  a  line  in  the  time- 
frequency  plane  [Wood  and  Barry,  1994]. 

In  this  paper,  we  investigate  the  application  of  the  WHT 
or  WRT  and  their  variants  to  GPS  post-correlation 
integration.  For  an  acquisition  scheme  wherein  frequency 
search  step  is  500  Hz  and  code  search  step  is  a  half  chip, 
GPS  despreading  correlation  is  typically  carried  out  over 
1  millisecond  (ms)  [Parkinson  and  Spilker  Jr.,  1996;  Tsui, 
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2000;  Misra  and  Enge,  2001;  Kaplan  and  Hegarty,  2006]. 
Such  integration  at  the  rate  of  1  kHz  can  tolerate  ±250  Hz 
errors  in  frequency  and  ±14  chips  in  code  phase  in  the 
sense  of  3  dB  loss  in  correlation  power.  In  many 
applications,  it  is  desired  to  further  integrate  the  1  ms 
correlations  over  a  rather  long  period  of  time.  When  the 
integration  interval  extends  to  and  beyond  20  ms,  the 
signs  of  data  bits  have  to  be  considered  in  addition  to 
changes  in  signal  frequency  [Yang  and  Han,  2006]. 

However,  the  combination  of  the  WVD  and  HT,  called 
the  Wigner-Hough  transform  (WHT),  or  the  combination 
of  the  WVD  and  RT,  called  the  Wigner-Radon  transform 
(WHT),  cannot  be  applied  directly  to  the  despread 
correlations  because  of  the  presence  of  unknown  data  bits. 
We  set  forth  a  block-implemented  WHT/WRT  using  the 
intra-block  conjugate  products  in  this  paper.  They  can  be 
used  to  jointly  estimate  the  initial  frequency,  chirping  rate 
(frequency  rate),  and  bit  sync.  After  introducing  a 
complex  model  for  post-correlation  GPS  signals,  we 
present  the  Wigner-Hough  transform  and  Wigner-Radon 
transform.  We  then  describe  a  variant  of  the  WHT  or 
WRT  based  on  intra-block  conjugate  products  together 
with  simulation  results  to  illustrate  the  operation  and 
performance.  Finally,  a  summary  is  provided  with 
concluding  remarks  and  future  work. 

COMPLEX  MODEL  FOR  POST-CORRELATION 
GPS  SIGNALS 

Consider  a  scheme  for  acquiring  GPS  signals  via  search 
in  time  and  frequency.  Assume  that  the  frequency  search 
step  is  500  Hz,  the  code  search  step  is  half  a  chip,  and  the 
integration  interval  of  the  despreading  correlation  is  1  ms. 
Such  integration  can  tolerate  ±250  Hz  errors  in  frequency 
and  ±14  chips.  The  resulting  complex  correlation  is  unique 
for  each  GPS  satellite  and  can  be  written  as: 

xn  =  b„A„  exp{j[27r(f0nTs  +  an2T2)+  <j>0 ]}  +  wn  (1) 

where  x„  is  the  complex  post-correlation  signal  with  the 
sample  index  n,  A„  is  the  amplitude,  is  the  initial  phase, 
Ts  =  1  ms  is  the  despreading  integration  interval  with  the 
sampling  rate  of  1  kHz,  b„  —  ±  1  is  the  unknown  data  bit ,f0 
is  the  initial  frequency,  a  is  the  chirping  rate  (sometimes 
a/2),  and  vv„  is  an  uncorrelated  noise  with  zero  mean  and 
unity  variance. 

For  GPS  C/A-codes,  the  baud  rate  is  50  Hz  and  the  data 
bit  thus  has  a  periodicity  of  20  ms.  There  will  be  20 
samples  per  data  bit  of  the  same  sign,  but  it  may  change 
sign  from  one  data  bit  to  the  next.  Knowing  which  sample 
is  the  start  of  a  data  bit  is  the  bit  synchronization  or  bit 
sync  process.  Without  knowing  the  bit  sync,  a  perfect 
correction  may  be  destroyed  if  there  is  a  data  bit  sign 
reversal  in  the  middle  of  integration.  It  is  equally 
important  to  determine  the  sign  of  data  bits,  from  which 
the  navigation  messages  can  be  demodulated. 


The  instantaneous  frequency  in  Eq.  (1)  represents  the 
frequency  error  between  the  incoming  signal  and  the  local 
carrier  replica,  which  can  be  written  as: 

Af„=f0nTs+cm2T2  (2) 


For  the  local  carrier  replica  fixed  at  a  search  frequency, 
the  complex  model  in  Eq.  (1)  is  useful  as  long  as  the 
frequency  error,  albeit  being  time-varying,  remains  within 
±250  Hz.  Otherwise,  the  search  switches  to  another  search 
frequency.  Over  a  time  interval  of  Th  the  maximum  chirp 
rate  for  a  fixed  carrier  replica  is: 


a  =  -sign(f0) 


500- 1 /o  I 

T, 


(3) 


where  sign(m)  is  a  sign  function,  so  that  the  change  in 
frequency  is  less  than  500  Hz  in  one  second.  The 
maximum  chirp  rate  is  therefore  a  =  500  Hz/s.  At  LI  = 
1575.42  MHz,  the  frequency  change  over  one  second  is 
related  to  the  underlying  acceleration  (in  g  =  9.8  m/s2)  by 
the  factor  of  51.49  Hz/g/s. 

The  signal  amplitude  An  includes  a  factor  R( A  r)  related  to 
timing  error  At,  where  /?(■)  is  the  correlation  function  for 
the  GPS  spreading  code,  and  a  factor  sin(nTsAf)  /  jiT  A/ 

related  to  frequency  error  A f  [Parkinson  and  Spilker, 
1996;  Tsui,  2000;  Misra  and  Enge,  2001;  Kaplan  and 
Hegarty,  2006].  As  a  result,  the  signal  amplitude  A„  is 
also  subject  to  variations  due  to  changes  in  A f.  Since  the 
change  is  rather  small,  the  amplitude  is  assumed  to  be 
constant  for  simplicity  in  this  paper. 

In  many  applications,  it  is  desired  to  further  integrate  the 
1  ms  correlations,  x„,  over  a  rather  long  period  of  time,  Ts. 
Because  of  the  unknown  data  bits  bn  and  varying 
frequency  errors  A f  it  is  counter-productive  to  sum 
complex  correlations  xn  directly.  Squaring  can  move  both 
the  data  bit  and  the  unknown  frequency  error  and  allows 
for  direct  summation.  However,  this  non-coherent 
integration  also  squares  noise  [Yang  et  al.,  2007]  and 
further  strips  off  all  valuable  information  about  the  signal 
that  is  needed  for  subsequent  processing  in  a  GPS 
receiver. 


WIGNER-HOUGH  AND  WIGNER-RADON 
TRANSFORMS 

The  post-correlation  GPS  signal  as  modeled  in  Eq.  (1)  is  a 
chirp  signal  with  a  linear  frequency  modulation.  The 
instantaneous  frequency  is  characterized  by  two 
parameters,  namely,  the  initial  frequency  f0  and  the 
chirping  rate  a.  The  Wigner-Ville  distribution  is  a  popular 
method  to  represent  the  signal  energy  in  the  time- 
frequency  plane  [Flandrin,  1999].  It  obtains  the  signal 
energy  distribution  by  applying  the  Fourier  transform  to 
the  centered  autocorrelation  function  (CAF)  of  a  time- 
varying  signal. 

For  a  continuous-time  signal  x(t),  its  Wigner-Ville 
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distribution  or  WVD  is  computed  as: 


Lxy\  (x,  y)  =  (-00  <  x  <  co,  y  =  sin  1 0(7- -  xcos0))  (7) 


W(t,f)  =  \x{t  +  -)x\t--)e~il,fTdT  (4) 

—oo  “  “ 

where  *  stands  for  complex  conjugate.  A  window 
function  is  typically  applied  so  as  to  restrict  the  integral 
limits  to  a  finite  length.  With  t  =  nTs,  z=  mTs,  and  /  = 
k/(N+l)Ts,  the  discrete  version  of  the  Wigner-Ville 
distribution  (DWVD)  can  be  written  as: 

N/2 

W(n,k)  =  2  2>(«  +  ni)x' (n  -  >n)exp(-j4mnk  /(N  + 1))  (5) 

m=—N  /  2 

where  the  total  number  of  samples  analyzed  is  /V+1 . 

In  addition  to  discrete  time  steps,  the  discrete  distribution 
also  has  discrete  frequency  steps.  This  may  be  perceived 
as  a  limitation  when  there  are  closely  spaced  signals  to 
resolve.  For  a  fixed  sampling  rate,  one  way  to  reduce  the 
frequency  step  is  to  zero-pad  the  data  samples  beyond  N. 
Interpolation  can  also  be  used  to  obtain  an  “increased” 
sampling  rate.  This  is  a  tradeoff  between  computation  and 
resolution.  The  selection  of  a  window  function  is  another 
important  design  issue  that  can  be  used  to  suppress  cross 
terms  when  multiple  chirp  components  are  present  in  the 
signal. 

The  Hough  transform  (HT)  is  a  popular  digital  image 
processing  technique  widely  used  to  find  a  straight  line 
from  an  image  [Duda  and  Hart,  1972].  For  a  polar 
coordinate  system  with  its  origin  at  the  center  of  the 
image,  a  line  can  be  represented  by  two  parameters,  r>  0 
and  0  e  [0,  2ji),  where  0  is  the  angle  between  the  x-axis 
and  the  origin-passing  perpendicular  to  the  line  and  r  is 
the  length  from  the  origin  to  the  intercept  of  the 
perpendicular  to  the  line  (the  closest  point).  The  line  can 
be  written  as: 

Tfl-:  (0,  r)  =  (0,  7-  =  xcos0  +  ysin0)  (6) 

For  a  single  point  (x,  y )  in  the  image  plane,  there  are  an 
infinite  number  of  lines  go  through  this  point,  each 
corresponding  to  a  sinusoidal  curve  given  by  Eq.  (6)  in 
the  (0,  r)  plane.  For  a  set  of  points  that  form  a  straight 
line,  their  sinusoidal  curves  cross  at  the  parameters  for 
that  line  in  the  (0,  r)  plane.  The  number  of  curves  that 
cross  at  the  same  point  in  the  (0,  r)  plane  is  the  number  of 
the  image  points  that  lie  along  the  same  line.  Therefore, 
the  HT  thus  converts  a  line  in  the  original  image  plane 
into  a  point  in  the  parameter  plane,  while  integrating  the 
image  values  along  the  line. 

This  description  makes  the  Hough  transform  to  be 
conceptually  very  close  to  the  Radon  transform  (RT) 
[Deans,  1983].  In  general,  the  RT  calculates  the  integral 
of  a  function  f[x,  y)  over  the  set  of  all  lines,  each 
characterized  by  two  parameters  in  the  polar  coordinates 
(0,  /').  In  the  x-y  plane,  a  line  is  represented  by  the  set  of 
all  points  along  the  line: 


Although  the  integration  can  be  carried  in  the  x-y  plane 
directly  using  Eq.  (6),  it  is  not  convenient  particularly 
with  possible  singularities  (0  =  90  or  270  degrees).  An 
easy  implementation  consists  of  rotate  the  x-y  plane  by  0 
into  the  d-r]  plane  and  then  integrates  along  the  line, 
which  is  now  parallel  to  the  77-axis  defined  as  ( E,  =  r,  -co 
<  77  <  00),  as  illustrated  in  Figure  2.  The  Radon  transform 
can  now  be  written  as: 

Rf(9,r)=  |  f(x,y)dl 

00 

=  | /(£cos0-77sin0,£sin0  +  77cos0)(777’  £=r  (8) 

—00 

It  has  been  shown  that  the  Radon  transform  performs 
better  than  the  Hough  transform  to  detect  lines  when  an 
image  is  corrupted  by  speckle  because  noise  can  be 
effectively  attenuated  in  the  Radon  transform  with 
explicit  summation  along  lines. 

Clearly,  both  the  HT  and  RT  also  involve  selection  of 
discrete  steps  for  r  and  0.  There  will  be  quantization 
errors,  which  again  is  a  tradeoff  between  computation  and 
resolution. 


Fig.  2  -  Integration  Path  for  Hough  and  Radon  Transfomis 

When  the  Hough  transform/the  Radon  transform  is 
applied  to  the  Wigner-Ville  distribution,  the  Wigner- 
Hough  transform/the  Wigner-Radon  transform  results, 
which  can  integrate  the  energy  belonging  to  a  chirp  signal 
along  a  line  in  the  time-frequency  plane.  The  chirp 


parameters  are  therefore  obtained  as: 

(/o ,  a)  =  arg  max  [  W (t,  f0  +  at)dt 

fa, a  J 

(9a) 

=  || x(t  +  -^-)x*(f  -  ■^)e~J2’TU°+a,) dtd r 

(9b) 
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where  both  integrals  are  limited  to  some  observation 
interval. 

Since  the  WVD  is  a  density  function  of  signal  energy  over 
the  time-frequency  plane,  when  integrating  along  all  lines 
in  the  WD,  Eq.  (9)  provides  the  maximum  likelihood 
(ML)  estimates  of  the  linear  FM  signal  parameters  in 
white  Gaussian  noise. 

The  basic  ideas  of  the  Wigner-Ville  transform  and  the 
Flough  transform  or  the  Radon  transform  described  in  this 
section  are  applied  in  the  following  sections  for 
integrating  the  post-correlation  GPS  signal  subject  to 
unknown  Doppler  changes  and  data  bits. 

WHT/WRT  WITH  INTRA-BLOCK  CONJUGATE 
PRODUCTS 

As  modeled  in  Eq.  (1),  the  post-correlation  signal  is 
available  at  a  rate  of  1  kHz.  Our  goal  is  to  integrate  a 
large  number  of  post-correlation  signals  x„,  for  n  =  1,  ..., 
N,  in  the  presence  of  the  unknown  frequency  f0,  chirping 
rate  a,  and  data  bits  b„.  Consider  an  integration  interval  of 
1  s  when  N  =  1000  for  Ts=  1  ms.  To  deal  with  unknown 
data  bits,  the  N  consecutive  signal  samples  are  grouped 
into  K  blocks;  each  block  will  contain  M  complex 
correlations  such  that  N  =  KM.  A  block  of  data  is  denoted 
by  a  vector  ^  =  [x20(i_l)+„  i  =  1  ,...,M]  =  [xu,  i  =  1  . 

When  M  =  20,  a  block  will  last  20  ms,  which  is  the 
duration  of  a  data  bit  for  GPS  C/A-codes. 

To  facilitate  presentation,  we  omit  the  terms  related  to 
noise  in  the  derivation.  Referring  to  Fig.  3,  we  re-label  the 
time  index  of  each  sample  in  a  block  relative  to  the  block 
center  by  t  ±  z/2  where  t  is  the  time  of  the  block  center  in 
units  of  MTS  and  z/2  =  0.5,  1.5,  ...,  9.5  in  units  of  Ts.  The 
centered  autocorrelation  between  two  samples  with  delay 
ris  given  by: 

z(t,  z)  =  x(t +  ^)x*  (t —  ^)  (10a) 

_  ^2gj2>r(/o+2a«)r 

t  =  ±  1,  ±3,  •••,  ±M  — 1,  t  =  \,  2,  •••,  K  (10b) 


nth  Block  (over  1  Data  Bit) 

Ts 

i  = 

1  2  3  .  8  9  !  10  11  .  20 

'/2  ‘/2 

„  l'/2< - * - >VA  .  n 

t  -  t/2  :  t  -  x/2 

9/2 

T  9V2 

Fig.  3  -  Time  Diagram  for  Intra-Block  Conjugate  Products 

To  arrive  at  Eq.  (10b),  it  is  assumed  that  the  two  samples 
share  the  same  data  bit,  that  is,  b(t+T/2)b(t-Tl2)  =  (±1)2  = 


1 .  The  aspect  of  bit  sync  will  be  addressed  later  in  this 
section. 

If  the  signal  is  strong  enough,  the  fast  Fourier  transform 
(FFT)  can  be  applied  to  z(t,  z)  over  z  for  each  block  t  and 
the  results  are  denoted  by: 

Z(t,f)  = 

Jr{[(z(t,r),  r  =  -M  +  l,...,-3,-l,l,3,...,M  —  1),  0^]} 

(11) 

where  J^{"}  stands  for  the  FFT  applied  along  the 
dimension  z,  0„  is  a  vector  of  n  zeros,  and  /u  is  the  factor 
of  M  zeros  padded  to  z(t,  z)  prior  to  the  FFT. 

The  peak  location  of  Z(t,f)  will  provide  an  estimate  of  the 
quantity  /o  +  2  at.  Repeating  this  operation  for  all  blocks 
yields  the  time-frequency  distribution  of  the  signal 
energy,  which  forms  a  straight  line  starting  at  f0  and 
moving  at  the  slope  of  2 a  as  shown  in  Fig.  4.  This  is  in 
fact  the  basic  idea  behind  the  Wigner-Ville  transform 
given  in  Eq.  (4),  which  then  requires  search  in  two 
unknowns,  f0  and  a,  while  integrating  the  signal  energy 
using  either  the  Hough  transform  or  Radon  transform  as 
discussed  in  the  previous  section.  As  shown  in  Fig.  4,  the 
estimates  of  the  unknowns  are  taken  as  those  that 
maximize  the  resulting  sum  S(fo,  «): 

(/0,a)  =  argmax.ST/0,a) 

fa, a 

S(f0,a)=  Ez(?’/) 

/=/0+2etf 

reP,..,*] 

where  the  arguments  are  converted  from  r  and  0  to  f0  and 
a. 

Alternatively,  the  FFT  can  be  taken  on  z(f,  r)  over  t  for 
each  z,  denoted  by: 


(12a) 

(12b) 


xn=x(t  =  nTs),n  =  l,  ...,N 


M  2 M 


KM 


J] 


FFT  over  t 


l  =  lax 

FFT  over  Search  Path  for /0 


S(f0,  cc). 


\  a 


->fa 


Fig.  4  -  Relationship  between  Various  Transforms 
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Fig.  5  -  Twenty  Delayed  Sums  for  Bit  Sync 


Z(l,T)  =  T,{[{z(t,T),  t  =  1,2,..., K),  0^]}  (13) 

Alternatively,  the  FFT  can  be  taken  on  z(t,  t)  over  t  for 
each  r,  denoted  by: 

Z(l,r)  =  J,{[(z(t,T),  t  =  0^]}  (13) 


As  shown  in  Fig.  4,  the  peak  location  of  Z(l,  r)  within 
+( k+  1  )K!2  is  obtained  as: 

/(r)  =  argmax  |  Z(/,r)  |  (14) 


l(r)  provides  an  estimate  of  the  quantity  2 ax.  With  r 
known,  the  chirping  rate  a  can  be  estimated  directly  as: 


a(r)  = 


1  kr) 

2  tT?N(k  + 1) 


(15) 


Another  advantage  of  this  approach  is  that  the  number  of 
blocks  K  can  be  increased  as  more  data  become  available 
(and  for  possible  recursive  implementation)  whereas  the 
number  of  delays  M  is  limited  per  data  bit  duration.  When 
the  signal  is  strong  enough,  the  peak  values  Z(  \{j) ,  r)  = 
KA2eJ 2*r  can  be  further  FFT-analyzed  over  t  and  the 
location  of  the  peak  of  this  second  FFT  provides  an 
estimate  of  the  initial  frequency  fo. 

When  the  signal  is  not  strong  enough  to  ensure  reliable 
detection  of  peaks  of  in  Z(/,  zj,  search  for  the  signal  path 
(/,  t)  is  required.  Again,  there  are  quantization  errors  in 
defining  the  search  path,  which  may  be  reduced  to  some 
extent  by  zero-padding.  For  each  a.  the  peak  location  /  is 
a  well-defined  linear  function  of  r.  As  a  result,  the  search 
is  one  dimensional.  The  FFT  is  used  to  integrate  along  the 
selected  path,  leading  to  S(f0,  a)  given  by: 

S{f0,a)  =  Tt{[Z(1,t),  1  =  2cct,  t  =  1,.,.,/AT]}  (16) 


where  the  arguments  are  converted  to  fo  and  a  as  shown 
in  Fig.  4.  The  signal  parameters  are  estimated  using  Eq. 
(12a)  as  those  that  maximize  S(fo,  a )  defined  in  Eq.  (16) 

In  the  above  derivations,  it  is  assumed  that  the  two 
samples  of  a  block  involved  in  the  conjugate  product 
share  the  same  data  bit  sign  such  that  b(t+T/2)b(t-T/2)  = 
(±  1  )2  =  1.  This  implies  that  this  block  is  perfectly  aligned 
(or  bit  synchronized  -  bit  sync)  with  data  bits.  For 
practical  implementation,  however,  the  parameter 
estimation  has  to  be  performed  at  the  same  time  as  the  bit 
sync.  One  simple  approach  to  bit  sync  is  to  maintain 
twenty  sums  of  blocks.  The  start  of  blocks  of  a  sum 
corresponds  to  a  possible  bit  transition  location  within  20 
ms  and  the  block  in  one  sum  is  delayed  by  one  sample  (1 
ms)  relative  to  the  preceding  sum,  as  illustrated  in  Fig.  5. 

The  kth  block  in  the  b'h  sum  is  therefore  constructed  as: 

y_k  =  C-*-20(jt— i)+i+*-i  ’  * =  i>— >20]  > 

b=l,  ...,20,k=l,  ...,K  (17) 

The  method  of  Eq.  (16)  can  be  applied  to  each  of  the 
twenty  sums  over  K  blocks,  denoted  by  Sh(  vb,  n  =  1,  ..., 

K).  The  sum  that  produces  the  maximum  value  provides 
an  estimate  of  the  data  bit  transition  (i.e.,  the  tentative  bit 
sync)  as: 

b  =  argmaxSb (yb  ...,yb  ...,yb  )  (18) 

b  — 1  — «  —K 

SIMULATION  RESULTS  AND  ANALYSIS 

Simulation  examples  are  presented  below  to  illustrate  the 
operation  and  performance  of  the  method  using  intra¬ 
block  conjugate  products  for  joint  signal  parameter 
estimation  and  bit  sync.  In  the  first  simulation  example, 
we  consider  an  ideal  case  where  there  is  no  noise  and  no 
data  modulation.  The  integration  interval  is  1  s  with  N  = 
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1000  for  Ts  =  1  ms.  The  data  are  divided  into  K  =  50 
blocks,  each  with  M=  20  samples.  The  initial  phase  is  9= 
32°,  the  initial  frequency  is/)  =  -250  Hz,  and  the  chirping 
rate  is  a/2  =  250  Hz/s. 

Fig.  6  shows  the  amplitude  spectrum  of  the  FFT  applied 
to  the  intra-block  conjugate  products,  Z(I,  r),  as  given  in 
Eq.  (13)  where  k=  3.  The  length  of  the  FFT  is  (k+\)K  = 
4x50  =  200  but  the  peak  value  is  50,  equal  to  the  length  of 
non-zero  data  points  K  =  50.  There  are  two  groups  of 
curves  and  the  locations  of  these  curve’s  peak  correspond 
to  2arfor  r=  ±0.5,  ±1.5,  ...,  ±9.5.  Because  of  the  sign  of 
t,  the  two  groups  appear  symmetrically  in  positive 
frequency  bins  (left)  and  negative  frequency  bins  (right). 
For  this  ideal  case,  the  estimated  chirping  rate  using  Eq. 
(15)  is  ia(r)=  250,  correct  for  all  t,  as  shown  in  Fig.  7. 
To  save  computation,  however,  only  one  group  of  curves 
is  in  fact  needed,  either  with  positive  or  negative  r. 

The  FFT-analysis  of  the  peaks  selected  from  Z(/,  r)  for  r 
=  1,  ...,  juM  is  equivalent  to  integrating  along  the  signal 
path.  As  shown  in  Fig.  8,  the  peak  location  corresponds  to 
the  initial  frequency/).  For  this  case  with  ju=  3  and  M  = 
20,  the  FFT  length  is  (/r±l)M  =  80.  The  frequency 
estimate  corresponding  to  the  peak  at  41  is  (41-80)*6.25  = 
-243.75  Hz,  which  is  quite  close  to  the  true  value  off0  =  - 
250  Hz  without  interpolation. 

The  FFT  integration  in  Fig.  8  is  done  along  a  special  path, 
which  is  determined  by  picking  the  peaks  at  (/(r),  r),  r  = 
1,  ...,  /uM.  Without  peak  detection,  the  signal  path  can  be 
searched  using  Eq.  (16).  The  results  are  shown  in  Fig.  9. 
For  this  ideal  case,  the  estimated  chirp  rate  is  250.00  Hz/s 
(vs.  all  =  250  Hz/s),  the  initial  frequency  at  the  peak 
value  is  -237.50  (vs./)  =  -250  Hz),  and  the  peak  value  is 
947.92  (vs.  the  full  value  of  1000)  without  interpolation. 

In  the  second  simulation  example,  all  the  conditions  are 
kept  same  as  in  the  first  simulation  example  except 
complex  noise  is  added  such  that  the  resulting  signal  to 
noise  ratio  is  SNR  =  -3  dB.  Under  the  Nyquist  rate,  the 
signal  passes  through  a  low  pass  filter  with  an  equivalent 
bandwidth  of  ±fJ2.  For  Ts  =  0.001  s,  SNR  =  -3  dB  results 
from  A  =  0.7,  which  also  corresponds  to  C/N0  =  27  dB- 
Hz. 

Fig.  10  shows  the  amplitude  spectrum  of  the  FFT  applied 
to  the  intra-block  conjugate  products  along  k.  Compared 
to  the  ideal  case  in  Fig.  6,  there  are  only  few  peaks  that 
can  be  reliably  detected  for  this  low  SNR  case.  But  the 
spectrum  is  still  symmetric,  indicating  that  only  half  of 
delays  (either  positive  or  negative  r)  would  be  needed. 
The  use  of  more  delays  did  not  necessarily  reduce  noise. 

It  is  then  not  a  surprise  to  see  from  Fig.  11  that  the 
chirping  rate  estimates  per  r  are  quite  noisy  except  for 
those  strong  peaks  in  Fig.  10.  In  particular,  those 
estimates  for  small  f  s  are  erroneous.  Nevertheless,  Fig. 


12  shows  the  FFT  integration  along  the  path  determined 
by  the  peaks  of  Fig.  10,  which  produces  the  initial 
frequency  estimate  rather  correctly.  As  compared  to  Fig. 
8,  the  noise  floor  is  however  raised  and  the  peak  is 
lowered  in  Fig.  12. 

Instead  of  picking  peaks  from  Fig.  10,  which  is  not 
reliable  for  low  SNR,  search  is  done  over  possible 
chirping  rates,  each  specifying  a  particular  path  to 
integrate  using  Eq.  (16).  The  results  are  shown  in  Fig.  13. 
Compared  to  Fig.  9,  the  noise  floor  is  raised  and  the  peak 
is  lowered.  Without  interpolation,  the  estimated  chirp  rate 
is  237.00  Hz/s  (vs.  all  =  250  Hz/s),  the  initial  frequency 
at  the  peak  value  is  -237.50  (vs.  f0  =  -250  Hz),  and  the 
peak  value  is  394.41  (vs.  the  value  of  947.92  in  the  noise- 
free  case). 

In  the  third  simulation  example,  data  bits  are  modulated 
onto  the  signal  samples.  In  the  example,  the  first  bit 
transition  occurs  between  the  16th  and  17th  samples.  In 
other  words,  the  16th  sample  is  the  end  of  the  first 
incomplete  data  bit  and  the  17th  sample  is  the  start  of  a 
second  data  bit.  Fig.  14  shows  the  twenty  delayed  sums  at 
K- 1  for  the  ideal  case,  from  which  Eq.  (18)  can  be  applied 
to  determine  the  bit  sync  at  17th.  Fig.  15  shows  the  twenty 
delayed  sums  for  SNR  =  -3  dB.  The  presence  of  noise  has 
not  lowered  much  the  peaks  for  those  sums  that 
hypothesized  bit  transitions  near  b  =  7.  For  this  example, 
sums  starting  at  those  samples  placed  the  bit  transitions  in 
the  middle  of  integration,  thus  rending  it  noise  like.  This 
explains  why  adding  another  noise  term  did  not  change 
much.  However,  the  peaks  for  those  sums  that  are  near  b 
=  17  have  been  reduced  in  magnitude  nearly  by  half  for 
SNR  =  -3  dB.  Nevertheless,  it  is  clearly  possible  and 
reliable  to  detect  the  cumulated  peak  corresponding  to  the 
true  bit  transition,  thus  achieving  bit  sync. 

The  product  of  two  samples  within  a  block  effectively 
eliminates  the  sign  of  data  bits,  thus  making  integration 
across  blocks  possible.  However,  other  means  have  to  be 
used  to  determine  the  data  bit  signs  from  which  to 
demodulate  any  message  that  may  be  embedded.  One 
method  using  inter-block  conjugate  products  that  allows 
for  joint  signal  parameter  estimation  and  bit  sync  as  well 
as  bit  sign  determination  will  be  presented  in  [Yang  et  ah, 
2008], 

CONCLUSIONS 

In  this  paper,  the  GPS  signal  after  the  1  ms  despreading 
integration  was  modeled  as  a  modulated  chirp  signal.  This 
model  is  useful  even  in  the  acquisition  mode  where  search 
is  conducted  at  a  frequency  step  of  500  Hz  and  a  code  lag 
of  Zi  chips.  In  this  model,  the  binary  modulation 
represents  the  data  bits  of  navigation  message  while  the 
chiip  accounts  for  changes  in  residual  Doppler  frequency. 
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To  further  integrate  the  signal,  a  new  method  was 
presented  based  on  the  basic  ideas  of  the  Wigner-Ville 
transform  (delayed  conjugate  products)  and  the  Hough  or 
Radon  transform  (integration  along  a  linear  path).  The 
method  was  based  on  the  intra-block  conjugate  products 
and  so  designed  for  jointly  estimating  the  signal 
parameters  (i.e.,  the  initial  frequency  and  chirping  rate) 
and  data  bits  (bit  sync).  Simulation  results  were  presented 
and  illustrated  the  operation  and  performance  of  the 
method. 

The  block-integrating  methods  presented  in  this  paper  can 
handle  large  frequency  changes  as  caused  by  relative 
motion  and  clock  drift.  It  can  also  remove  data  bits  in  an 
open-ended  manner,  allowing  for  new  samples  to  be 
added  on  in  the  integration  process  as  they  become 
available.  These  and  other  features  make  it  suitable  for 
such  applications  as  in  high  dynamics  and/or  weak  signal 
environments.  Beyond  naturally  occurring  changes  in 


Doppler  frequency,  modulated  chirping  signal  may  be 
used  for  signaling  in  special  CDMA  applications  and 
MIMO  radar  for  detection  of  stealthy  targets  among 
others. 
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